Improving Noise Robustness in 
Subspace-based Joint Sparse Recovery 

Jong Min Kim, Ok Kyun Lee and Jong Chul Ye 



(N 

^ ■ Abstract 

In a multiple measurement vector problem (MMV), where multiple signals share a common 
sparse support and are sampled by a common sensing matrix, we can expect joint sparsity to 
enable a further reduction in the number of required measurements. While a diversity gain from 
. joint sparsity had been demonstrated earlier in the case of a convex relaxation method using an 

I1/I2 mixed norm penalty, only recently was it shown that similar diversity gain can be achieved by 
greedy algorithms if we combine greedy steps with a MUSIC-like subspace criterion. However, the 
main limitation of these hybrid algorithms is that they often require a large number of snapshots or 
a high signal-to-noise ratio (SNR) for an accurate subspace as well as partial support estimation. 
O ■ One of the main contributions of this work is to show that the noise robustness of these algorithms 

can be significantly improved by allowing sequential subspace estimation and support filtering, 
P\J ■ even when the number of snapshots is insufficient. Numerical simulations show that a novel 

^ . sequential compressive MUSIC (sequential CS-MUSIC) that combines the sequential subspace 

estimation and support filtering steps significantly outperforms the existing greedy algorithms and 
t:;;^- ■ is quite comparable with computationally expensive state-of-art algorithms. 
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I. Introduction 

We study a multiple measurement vector (MMV) problem, where multiple signals share a same 
common sparse support set and each signal is measured by multiplying it by a measurement matrix. 
An MMV problem is one way in which multiple correlated signals can appear in a signal ensemble, 
and MMV problems also have many important applications 0], 111, 111, lU. A central theme in 
these studies has been that joint sparsity within signal ensembles enables a further reduction in 
the number of required measurements Q, iQ, where the number of measurements required per 
sensor must account for the minimal features unique to that sensor Q, HI, 191, ifTOll . ifTTI . |[T2ll . 
Indeed, for the case of an h/h mixed norm approach, Obozinski et al. flT| showed that a near 
optimal diversity gain can be achieved. 

Recently, Kim et al. Q and Lee et al. |12| independently showed that such a diversity gain 
can be also achieved in a new class of greedy algorithms by exploiting the so-called generalized 
(or extended) multiple signal classification (MUSIC) criterion Q, |[T2l . More specifically, these 
algorithms obtain a partial support estimate using a conventional MMV greedy algorithm, and 
then the atoms corresponding to the partial supports are augmented into a data matrix to obtain 
an augumented signal subspace estimate. Finally, a MUSIC-like |[T4l criterion is derived for the 
augmented subspace to find the remaining support. The hybridization makes these hybrid greedy 
algorithms fully utilize a diversity gain so that the algorithms outperform all the existing greedy 
methods. 

The performance improvement of these greedy algorithms is substantial and nearly achieves the 
Iq bound when a signal subspace and partial support estimation are accurate due to a sufficient 
number of snapshots or high signal to noise ratio (SNR) Q, |[T2ll . However, if either of these 
estimation is erroneous owing to an insufficient number of snapshots or low SNR, performance 
degrades. Similar observations have been made in the literature on classical array signal processing 
ifTSl . lfT6l . In ifTSl . prior knowledge of the direction-of-arrival (DOA) has been incorporated to 
improve the performance of MUSIC by filtering out the known sources via orthogonal projections. 
However, as shown in [|16il . such orthogonal projection is suboptimal from a statistical standpoint. 

While increasing the number of snapshots is relatively easier in classical sensor array signal 
processing problems, in some MMV problems such as parallel MR imaging an additional 
snapshot requires a hardware change by adding a new receiver coil. Hence, in these problems, 
exploiting other dimensions would be beneficiall. We are aware that joint sparse recovery methods 
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such as Bayesian approaches ifTTl . ifTSl . or convex optimization techniques |fT9l . are shown to be 
statistically robust in the direction of arrival estimation problems, as first demonstrated by Malitov 
et al. II20I and further developed by Stoica et al. |18]. However, these approaches are usually 
computationally expensive for MMV problems with a large number of sensors, so we need a 
new greedy algorithm that achieves a similar optimal performance with a significantly reduced 
computational complexity. 

Therefore, one of the main goals of this paper is to address how these hybrid greedy methods can 
be made robust without increasing the number of snapshots. One important contribution is a new 
theory explaining that the generalized MUSIC criterion is a special case of a new subspace criterion 
that can be used to derive two sequential strategies to improve the accuracy of an augmented signal 
subspace estimation. More specifically, a forward greedy subspace estimation step improves the 
robustness of an augmented signal subspace estimation by adding newly discovered atoms in the 
MUSIC step, whereas the backward support filtering provides additional robustness by eliminating 
the inaccurate portion of support estimates. By combining the two steps, we develop a novel 
sequential CS-MUSIC algorithm that is robust, even with a limited number of snapshots. Using 
theoretical noise analysis as well as numerical simulation, we show that the sequential CS-MUSIC 
is superior to the existing subspace-based greedy algorithms and exhibits similar performance 
behavior to the mixed-norm |[T3]| . ll20l or Bayesian approaches ifTTll with a significantly reduced 
computationally complexity. 

II. Generalized Subspace Criterion 
A. Notations and Mathematical Preliminaries 

Throughout the paper, x* and correspond to the z-th row and j-th column of matrix X. When 
/ is an index set, , Aj corresponds to a submatrix collecting corresponding rows of X and 
columns of A, respectively. The rows (or columns) in M" are in general position if any n collection 
of rows (or columns) are linearly independent. 

Definition 1 (Canonical form noiseless MMV [7]): Suppose we are given a sensing matrix A G 
j^mxn ^j^^ observation matrix B G M™^*^ such that B = AX^ for some X* G M"^*^ and 
||X*|| = |suppX| = k, where m, n, and r are positive integers (r < m < n) that represent the 
number of sensor elements, an ambient space dimension, and the rank of an observation matrix, 
respectively. A canonical form noiseless multiple measurement vector (MMV) problem is given as 
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an estimation problem of A;-sparse vectors X G R"^'' using the following formula: 

minimize ||X||o (1) 
subject to B = AX, 

where ||X||o = |suppX|, suppX = {1 < -i < n : x* 7^ 0}, x* is the i-th row of X, and the 
observation matrix B is full rank, i.e. rank(i?) = r < k. 

Recall that every MMV problem can be converted to a canonical form MMV by using a singular 
value decomposition and dimension reduction as described in |7 |. Hence, in this section, we assume 
that an MMV problem assumes the canonical form. However, this assumption will be relaxed later 
in noise analysis. 

B. Generalized Subspace Criterion 

Note that the generalized MUSIC criterion in Q requires < 52f^_^_^i{A) < 1. This implies 
that, if a sensing matrix is obtained from a random Gaussian and if measurement is noiseless, then 
we have the following minimal sampling condition Q: 

m > (1 + 5){2k — r + 1) for some 6 > 0. 

If we have a redundant sampling m ^ 2k — r + 1, the following theorem can be used instead as 
the extension of the generalized MUSIC criterion in Q. 

Theorem 1: Suppose 1 < / < r and we have a canonical MMV model AX = B with a sensing 
matrix A that satisfies an RIP condition with < 52fc-r+i(^) ^ ^- Furthermore, suppose the 
nonzero rows of X are in general position. Then, for a given index set / C { 1 , • • • , n} such that 
|/| < min(2(/c — r) + /, k) and \I \ suppXj < k — r + I, the following statements are equivalent: 

(i) |/n suppX| >k-r + l] 

(ii) rank[^/ 5] < |/| +r. (2) 

Proof: (i)<;=^(ii): Assume that \I n suppXj > k — r + 1. Then \I\ + r < 2k — r + I < m 
since 62}^_^_^i{A) < 1. If we take / C (/ n suppX) such that \I\ = k — r + 1, then 

dim{R[Aj B]) < dim(yl5) = k. 

However, \I\ + r = k + I so that [Aj B] is not of full column rank. Hence [Aj B] is not also of 
full column rank since I C I. 
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Conversely, if we assume (ID, there are p G M'^' and q G M'" such that Ajp + AXq^ = and 
[P^q]^ / 0- If we let p e M" by p-^ = p and p^" = 0, we have A[p + Xq\ = 0. Since 
||p + Xq||o < |/\suppX| + |suppX| < 2A; — r + /, by the RIP condition, we have p + Xq = so 
that supp(p) = supp(Xq) C suppX. If we assume that \I n suppX| < k — r, then ||p||o < k — r 
but ll^qllo > A; — r + 1 since the nonzero rows of X are in general position. This is impossible 
so that |/n suppX| > A; - r + 1. ■ 

Note that the conditions |/| < min(2(/i;— r)+/, k) and |/\suppX| < k—r+l in Theorem[T]do not 
imply that there is a unique index set /; rather. Theorem [T] says that multiple index sets / can exist 
for a given /. For example, if / = 1, any index set / such that |/| = k—r+l, • • • , min(2(A;— r) + l, k) 
that satisfies the condition |/\suppX| < k — r + l, can be used to test conditions (i)-(ii) in Theorem 
[T] Furthermore, if we choose \I\ = k — r + 1, Theorem [T] is reduced to the following generalized 
MUSIC criterion in Q. 

Corollary 1 (Generalized MUSIC Criterion i[Zl/j-' Suppose we have a canonical MMV model 
AX = B with a sensing matrix A that satisfies an RIP condition with < SJ^k-r+i 

{A) < 1. 

Furthermore, suppose the nonzero rows of X are in general position. Then, for C suppX 
with \Ik-r\ = k — r and any j G {1, • • • ,n}\ Ik-r, we have j G suppX if and only if 

rank[Aj^_u{i} B] < k + 1 (3) 

or equivalently 

Proof: For an index set / such that |/| = k — r+l, the condition |/\suppX| < k — r+l always 
holds. Therefore, |/n suppX| = k — r + 1 is equivalent to / C suppX. Hence, for / = Ik-r U {j} 
such that Ik~r C suppX and \Ik-r\ = k — r, Eq. (|2]) is equivalent to Eq. which is equivalent 
to say Sij G R{[Ai^_^B]) or a*P^^j^^ b])^J ~ "^^^^ concludes the proof. ■ 

Remark I: If r = k, the conventional MUSIC criterion can be trivially derived. 

Remark 2: The subspace R{[Aj^_^ B]) is called augmented signal subspace. This name was 
first coined in | fT2| . 

So far, we have shown that Theorem [U can reproduce the existing results. However, one of the 
important byproducts of the theorem is the following form, which will be used extensively in the 
following sections. 
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Corollary 2: Suppose 1 < I < r and we have a canonical MMV model AX = B with a 
sensing matrix A that satisfies an RIP condition with < ^2k-r+i^'^^ ^ ^- Furthermore, suppose 
the nonzero rows of X are in general position. Then, for an index set / C { 1 , • • • , n} such that 

|/| < min(2(/c — r) + 1, k) and \I \ suppX| < /c — r + if we have |/ n suppXj = k — r + q iox 
some (/ > 0, then 



Proof: Take an Ik-r C (In suppX) with \Ik-r\ = k — r and let Jk-r = (in suppX) \ Ik-r, 
where |Jfc-r| = 1- Then by Theorem [T] we have 

Twk[Ai\j^_^ B] = \I\ Jk-r\ + r = \I\+r -q. 

Then for any j G Jk-r, G R{[Ai^_^^ B]) = R{As) so that 

rank[yl/ B] = Ta.-Dk[Ai\^j^_^ B\ = \I\ + r — q 



III. Sequential Compressive MUSIC Algorithm 

By employing the results in the previous section, this section first develops forward or backward 
greedy steps. Then, by combining the two approaches, we can derive a novel sequential CS-MUSIC 
algorithm. 

A. Forward Greedy: Sequential Subspace Estimation 

In 111, the CS-MUSIC first determines k — r indices of suppX with CS-based algorithms such 
as 2-thresholding or S-OMP, and then it recovers the remaining r indices of suppX using the 
generalized MUSIC criterion. For this, a projection operator onto the noise subspace is calculated 
as the orthogonal complement of the augmented signal subspace R{[Ai^_^ B]). However, the 
following result can further extend the existing generalized MUSIC criterion Q . 

Theorem 2: Suppose 1 < Z < r and we have a canonical MMV model AX = B with a sensing 
matrix A that satisfies an RIP condition with < 52i^_j.^i{A) < 1. Furthermore, suppose the 
nonzero rows of X are in general position. Then, if we have an index set / C {1, • • • , n} such 
that |/| < min(2(A; - r) + / - 1, A; - 1), |/ \ suppX| < A; - r + / - 1 and |/ n suppX| > k - r, 
we have for j ^ I, j ^ suppX if and only if 



rank[A/ B] = \I\ + r - q. 



(4) 



since R{[Aj B]) = R{[Ai\j^_^ B]). 



vank[Aj B] = rank[Aju{j} B] 



(5) 
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or equivalently 

Proof: By the condition we have \I U {j}\ < min(2(A; — r) + /, k) and | (/ U {j}) \ suppXj < 
k — r + I so that we can apply Corollary |2] for / and / U {j} since \I n suppX| > k — r. If 
\I n suppX| = k — r + q for some q > 0, then we have rank[A/ B] = \I\ -\- r — q. Then for any 
j ^ /, if we have j G suppX, |(/ U {j}) n suppX| = k — r + {q + 1) so that we have 

rankfA/ujj} = |/| + 1 + r - (q' + 1) = |/| + r - = rank[yl/ B]. 

On the other hand, if we have j ^ suppX, |(/ U {j}) H suppX| = k — r + q so that we have 

rank[A/u{j}] = |/| + l + r- g'> rank[^/ B]. 

Finally, dD is equivalent to aj G R{[Aj B]), which is also equivalent to This completes the 
proof. ■ 

Remark 3: Note that R{[Ai B]) = R{[Ai^_^ B]) and dimR{[Ai^_^ B]) = k for all / C suppX 
and k — r < \I\ < min {2{k — r) + 1 — 1, k — 1). This implies that we first need to find I^-r 
support using a compressive sensing algorithm, then we augment newly added supports into the 
initial estimate Ik~r- As will be shown later in noise analysis, such a greedy procedure improves 
the accuracy of the augmented signal subspace estimation. 

Remark 4: The greedy procedure can even be performed in a critically sampled case, i.e. I = 1. 
In this case, we can augment atoms up to min(2(A; — r),k — 1), which is always bigger than 
adding only k — r atoms. However, the number of possible augmentation increases with a redundant 
sampling, which makes the algorithm more robust. 

Theorem |2] leads us to the following sequential algorithm (SeqSubspace), as in Table I. Note 
that the algorithm can be combined with any joint sparse recovery algorithm that provides & k — r 
initial support estimate. 

B. Backward Greedy: Support Filtering 

As discussed before, we can easily expect that the performance of the generalized MUSIC 
step is highly dependent on the selection of A; — r correct indices of the support of X. Note 
that this is a very stringent condition. In practice, even though the first consecutive steps of, for 
example, S-OMP, may not provide aU true partial supports, it is more likely that among a A;-sparse 
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TABLE I 

Algorithm: /= SeqSubspace(A, B, Jfe_,.) 

- Set q = and / = Ik-r- 

- While g < r, do the following procedure: 

1. Perform an SVD of [Ai B] = [(7i, (7o]diag[Ei, SolfFi, I/q]*, 
where Ei = diag[CTi, ■ ■ ■ ,crk] and Eo = diag[(jfc+i, • ■ • ,ak+q] 
and (Ti > (72 > • • • > ak+q. 

2. Take = argmiiij^j \\PR(Ui)^jf- 

3. Set I :— I U {jq}, let g := g + 1 and goto step 1. 

- Return /. 



support estimate of S-OMP, part of the supports (not in sequential order) can be correct. In fact, 
an information theoretical analysis of a partial support recovery condition in single measurement 
vector CS (SMV-CS) ||2TI showed that the required SNR condition of a partial support recovery 
is much more relaxed than that for a full support recovery. Hence, if the estimate of the support 
of X has at least k — r indices of the support of X and we can identify them, then we can expect 
that the performance of the compressive MUSIC will be improved. When {f.^^) is small, we may 
apply the exhaustive search, but if both k — r and r are not small, then the exhaustive search is 
hard to apply so that we have to find some alternative method to identify correct indices from an 
estimate of suppX. 

Indeed, our new algorithm requires that k — r + 1 supports (not in sequential order) out of a larger 
support estimate is correct. Then, the location of a correct k — r support can be readily estimated 
using the following backward support filtering. Compared to a forward greedy procedure that 
improves the accuracy of the signal subspace estimation, the backward support filtering criterion 
can improve the accuracy of a partial support recovery, and, hence, the corresponding accuracy of 
an augmented signal subspace. 

Theorem 3 (Backward support filtering criterion): Suppose 1 < I < r and we have a canonical 
MMV model AX = B with a sensing matrix A that satisfies an RIP condition with < 
'^2ik-r+«(^) ^ ^- Furthermore, suppose the nonzero rows of X are in general position. Then, if we 
have an index set / C {1, • • • , n} such that |/| < min(2(A; — r) + Z, k), \I \ suppXj < k — r + I 
and 1 1 n suppX| > k — r + 1, then we have for j £ I, j £ suppX if and only if 

rank[^/\|j} B] = rank[Aj B], 
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or equivalently 

Proof: Assume that \I n suppXj = k — r + q, where q > I. Then, by Corollary |2l we have 
rank[A/ B] = \I\ + r — q. Noting that / \ {j} satisfies the assumptions of Corollary |2]for any 
j G /, if we have j G suppX, \I n suppXj = k — r -\- {q — I) so that we have 

rank[^7\{j} B] = \I\ - 1 + r - {q - 1) = \I\ + r - q = rank[^/ B]. 

On the other hand, if we have j ^ suppX, \I n suppXj = k — r + q so that we have 

Ta.nk[Aj\[jy B] = \I\-l + r- q=\I\+r-q-l< rank[A/ B]. 

Finally, due to the rank condition, we know slj G B]) if and only if j G suppX. Hence 

} B])^j = if and only if j G suppX. That completes the proof. ■ 
Theorem |3] informs us that if we have a partial estimate of support of X that has at least k — r + 1 
correct indices of support of X, we can identify the correct part of the estimated partial support 
of X by using the backward support filtering criterion as described in Table II. 

TABLE II 

Algorithm: Ik-r= SupportFiltering(y4,i3, /) 

- For all j £ I, calculate the quantities C^{j) := j 

- Making an ascending ordering of (^(j) for j £ /, choose indices that 

correspond to the first k — r indices and put these indices into Ik-r- 

- Return Ik-r- 



Remark 5: Due to the condition |/| < min(2(A; — r) + l,k) in Theorem [3l we can include k- 
sparse support estimate I in a support filtering step if r < (A; + l)/2. Note that this is always true 
regardless of r if / = A; or (^2^ < 1. However, if r > {k + l)/2, we can use the following heuristics. 
First, just include the first 2{k — r) + l support estimate of / for a support filtering. Since, in most 
greedy algorithms, the earlier greedy steps are more likely to succeed, correct k — r + I supports 
are more likely to be included. Hence, we can filter out the remaining indices j £ I such that 
j ^ suppX. 
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C. Sequential CS-MUSIC 

By combining the forward and the backward greedy steps, this paper develops the following 
sequential CS-MUSIC algorithm decribed in Table m. Note that this algorithm assumes that the 
sparsity level k is given as a priori knowledge. (The estimation problem of an unknown k will be 
discussed later.) 

TABLE m 

Algorithm: h= SeqCSMUSIC(^,y,fc,r) 
Input: k,r,A€ R""""" , Y e W^""^ 
Oulpul: /,■ support cslimalc i;, 

- Estimate the k support estimate Ik of suppX using any MMV algorithm. 

- U := Rank-r signal subspace estimate of R{Y). 

- Ik-r :=SupportFiltering(A,?7,7fc). 

- Ik :=SeqSubspace(yl, U, Ik-r)- 

- Return Ik- 



IV. Noisy Performance Analysis of Sequential CS-MUSIC 

A. Improving Noise Robustness Using Sequential Subspace Estimation 

In practice, measurements are noisy, so the theory we have derived for noiseless measurements 
should be modified. Suppose a noisy MMV model is given by: 

Y = AX + W, 

where Y G W^^^ are noisy measurements corrupted by an additive noise W G W^^^ , and 
N denotes the number of snapshots. Then, using singular value decomposition, we can find the 
following canonical MMV problem: 

S = AX + W, 

where S G M"*^^ is the rank-r signal subspace estimate of Y and r < N denotes the numerical 
rank of Y. Due to the noise, S is peturbed from the noiseless signal subspace S such that 
R{S) = R (y^^^ ' which leads to errors in the augmented signal subspace. The following theorem 
characterizes how much perturbation in an augmented signal subspace can be endured by a 
generalized MUSIC step. 
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Theorem 4: For < / < r, if we have Ik-r+i C suppX such that \Ik^r+i\ = k — r + I and 
singular value decomposition of [A/^^^ 5*] and [Ai^_^^^ S] as 

S] = C/iSi^i, [Ai,_^^, S] = ilitiV^ + UotoVo*, 

where Si = diag[ai, • • • ,ak], So = diag[afc+i, • • • ,ak+i] and cti > (72 > • • • > &k+i, then, for 
any j ^ suppX and q G suppX we have 

II a -11^ ^ II a 11^ 

provided that 

WPR^-PRiuJ <l-7. (7) 

and m,n — )• oo and 7 = liuin^ook/m < 1 . In other words, a generalized MUSIC step finds 
correct supports if Eq. Q is satisfied. 

Proof: Noting that .jajlp = for j G suppX by the generalized MUSIC criterion, for 

any j ^ suppX and q G suppX, we have 



||p± ||2 _\\p± „ ||2 



- ^*j^R{u,)^j ~ ll^m) ~ PR{u,)\\\\^q\? ■ (8) 

Since Ojj's are i.i.d. normal distribution with zero mean and variance 1/m and is independent 
of P'^fjj ^ for any j ^ suppX, sl^P^^^ ^slj is a chi-squared random variable of degree of freedom 
m — k since rank(C/i) = k. Also, for each 1 < j < n, m||aj|p is a chi-squared random variable 
with degree of freedom m so that we have, by Lemma 3 in [22], lim max HaJp/m = 1 since 

n—s-oo l<j<n 

lim„_j.oo (log n)/m = 0. Since lim„^oo(log {n — k))/{m — k) = 0, by Lemma 3 in ll22l . we have 
lim min ma*,Pj-,- ,a.j/(m — k) = 1 so that lim min a*:P:^,~ ,ao/ max llaJP = 1 — 7. 
Hence, Eq. ^ is positive provided that Eq. (|7]l holds in the large system limit. This completes the 
proof. ■ 
Therefore, by minimizing the perturbation in the augmented signal subspace \\Pj^0^-^ ~ Pr{Ui)\\^ 
we can make the generalized MUSIC step more robust. Unfortunately, the direct minimization of 
the perturbation of the subspace is not easy. Instead, we are interested in minimizing the following 
upper-bound of the perturbation, whose proof can be found in Appendix A: 
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where A = US' — 5*11 and crk{[Aj^_^^, S]) denotes the fc-th largest singular value of [Ai^_^^^ S]. 
Then, we have the following theorem: 

Theorem 5: Let / > 0. For I^-r+i C suppX such that = k — r + I, the generalized 

MUSIC steps find the remaining r — / support provided that 



CTfc([^/,_„„ S]) 



> 1 + 



1 



(10) 



A ' ' ' 1-7 

and m, n — )• oo and 7 = lim„_j.oo k/m < 1. 

Proof: This can be trivially proven by pugging Eq. Q in the inequality ([7]). ■ 
Note that for Ik-r+i C suppX, we have rank([A4^^, S]) = A; so that the set of columns 
of [Ai^_^,^^ S] is a frame in ii(^suppx) with lower frame bound (T^([Ajfc_r+i S\). In this case, 
as shown in Fig. [TJ ak{[Ai^^_^^^^ B]) is an increasing function of /, so as / increases, the frame 
becomes more redundant and the lower frame bound become larger. Hence, the left side of Eq. ([TOl l 
becomes larger. 
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I 



Fig. 1. (Jk [[Ak-r+i, S\) values with increasing I. Simulation parameters are n = 128, = 8, r = 6. 



This observation provides us an important error correction scheme. Note that the SNR condition 
Eq. (ITOl) is still the same even if we find a support index ji in a greedy manner as follows: 

= arg mill \\P^-.a.jf. (11) 

However, as SNR condition Eq. (fTOl ) is a sufficient condition, a non-zero probability of ji being 
in the true support exists even though Eq. (ITOl) is not satisfied. (This is especially true if we select 
only one index rather than choosing all r — i indices). If a correctly found index ji is augmented for 
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the next step of sequential subspace estimation, then it is more likely that the condition in Eq. ([TOb 
can be satisfied in the following greedy steps since the left side term of Eq. (fTOl) is an increasing 
function of I thanks to the inclusion of a correct index ji. As soon as an SNR condition is satisfied, 
the remaining greedy steps will succeed since the condition is sufficient. Therefore, even when 
a sufficient SNR condition is not satisfied initially, the proposed sequential subspace estimation 
technique exploits the possibility of finding a correct index to improve the noise robustness, which 
was not possible in an original MUSIC step. 

B. Improving Noise Robustness Using Support Filtering 

Using similar techniques, we can derive the following sufficient condition for the success of 
support filtering. 

Theorem 6: Let m > k + r. Suppose we have an index set / C {1, • • • , n} such that |/| = k, 
\Ic\ > k — r + I, where Ic ■= (suppX) n I. Then, we have 

mill ll-Ppr. eiajlP > max \\P^,a q^^j\\'^ (12) 

provided that 

^k{qk-r,I) 1 

A ^'+1-^(1 + a)' ^ ^ 

and m,n — )- oo, where a{q,I) = m.ax{ak{[Aj^ S]) : T C Ic \ {q}, 1^1 = ^ — for q e Ic 

and Ic = {qi,--- jQ\ic\}y which satisfies a{qi,I) > a{q2,I) > ••• > iT(g|/p|,/), and 7 := 

limn->oo fc/m, a := limn^oor/k. 

Proof: See Appendix B. ■ 

In Theorem [6l cJkiqk-r^ I) increases when an initial support estimation has more correct support 
because of the equation in which dk{q, I) is given by the maximum value out of possibilities. 
Moreover, if we increase the ratio m/{k + r), then the right-hand side of ([T3] ) decreases so that we 
can expect a greater possibility of accurate support filtering with an increased redundant number 
of samples than the critical sampling rate. 

To confirm a support filtering useful for performance improvement, we examine the cases 
where the sufficient condition for an initial support estimation is less favorable than that of a 
support filtering. Characterization of such cases should be done with respect to a particular initial 
support estimation algorithm. For example, in the case of a subspace S-OMP for an initial support 
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estimation, a sufficient condition for the success of subspace S-OMP for q > is given by iTTl i 



< 



a - a^{2 - F{a)) 



(14) 



where F{a) is an increasing function such that F{1) = 1 and limQ^o+ = 0, which is defined 

a/o'*^''"'* xdXi{x), dXi{x) = (-y/ (4 — x)x)/{2TTx)dx is the probability measure with 



as F(a) 



support [0,4], < ti{a) < 1 satisfies f^^^^"'^ dsi{x) = a, and dsi{x) = {1/tt)V4 — x'^dx is a 



probability measure with support [0,2]. Now, the SNR condition for support filtering in Eq. (1131) 
can be translated into a threshold of the allowable augmented subspace perturbation of 1 — 7(1 + 0;), 
when m > r + k. Hence, the gap between the two bounds is given by 

a - a^{2 - F{a)) 



/(7>a) 



(1-7(1 + a)). 



Fig. |2] characterizes the function /(7, a). In region A, f{'y,a) < and m > r + k, hence, the 
support filtering has more noise robustness and can correct errors from subspace S-OMP. As we can 
see from Fig. |2j support filtering is effective in most of the practical sampling rate, and especially 
when we have redundant samples or rank(y) is relatively small. The larger size of region A that 
favors support filtering again confirms that support filtering is a quite useful technique to improve 
noise robustness. 




Fig. 2. Feasibility region where support filtering becomes more effective than subspace S-OMP. 
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V. NUMERICAL RESULTS 

In this section, we perform extensive numerical experiments to validate the proposed algorithm 
under various experimental conditions, and compare it with respect to existing joint sparse recovery 
algorithms. 

A. Dependency on Snapshot Number 

First, we demonstrate that a sequential CS-MUSIC is less sensitive to the number of snapshots. 
The simulation parameters were as follows: m G {1, 2, . . . , 30}, n = 128, = 8, r = 4, and 
N G {6, 16, 256}, respectively. The elements of a sensing matrix A were generated from a Gaussian 
distribution having zero mean and variance of 1/m, and then each column of A was normalized to 
have an unit norm. An unknown signal X with rank(X) = r < k was generated using the same 
procedure as in |[T2l . Specifically, we randomly generated a support /, and then the corresponding 
nonzero signal components were obtained by 

= ^A^ , (15) 

where ^ G R^^*" and A were set to random orthonormal columns and the identity matrix, 
respectively, and <I> G W^^ were made using Gaussian random distribution with zero mean and 
variance of 1/A^. After generating noiseless data, we added zero mean white Gaussian noise to 
have SNR = 30dB measurements. We declared success if an estimated support was the same as 
a true suppX, and success rates were averaged over 1000 experiments. 

Fig. [3] shows success rates of a sequential CS-MUSIC compared to that of CS-MUSIC or SA- 
MUSIC. Since SA-MUSIC in |12| is equivalent to CS-MUSIC for a normaUzed A matrix, the 
original code of SA-MUSIC was used for fair comparison. As shown in Fig [3l sequential CS- 
MUSIC exhibits nearly similar recovery performance for various snapshot numbers, whereas the 
original form of CS-MUSIC/SA-MUSIC requires a large number of snapshots to achieve maximum 
performance. 

In order to identify the contribution of the forward and backward greedy steps in the performance 
improvement, we perform additional experiments using the same simulation setup. Fig. |4] illustrates 
the performances of sequential CS-MUSIC, a variation of sequential CS-MUSIC without backward 
support filtering, and the original CS-MUSIC/SA-MUSIC algorithm for = 6, 16, respectively. 
Here, an initial k - r support for CS-MUSIC/SA-MUSIC and the sequential CS-MUSIC were 
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Fig. 3. Snapshot dependent performance behaviour of the sequential CS-MUSIC and the original CS-MUSIC/SA- 
MUSIC. The simulation parameters are n = 128, fc = 8, r = 4 and SNR = 30dB. 



estimated using an identical subspace S-OMP algorithm in Q, |[T2l so that performance differences 
came only from the sequential subspace estimation step. For a bigger N where the signal sub- 
space error is small, performance improvement due to the sequential subspace estimation was not 
remarkable. However, the advantages of sequential subspace estimation is especially noticeable for 
a small number of snapshots where a subspace estimation is prone to error. On the other hand, the 
backward support filtering is beneficial for all ranges of snapshots since it corrects the contribution 
of a partial support estimation error in a subspace S-OMP step. 

B. Performance Comparison with State-of-Art Joint Sparse Recovery Algorithms 

To compare the proposed algorithm with various state-of-art joint sparse recovery methods, the 
recovery rates of various state-of-art joint sparse recovery algorithms such as CS-MUSIC/SA- 
MUSIC, /1//2 mixed norm approaches lEl, EQI, ES, and M-SBL 113, are plotted in Fig. |5] 
along with those of a sequential CS-MUSIC. Among the various implementation of mixed norm 
approaches, we used high performance SGPLl software |[23l . which can be downloaded from 



http : // www.cs.ubc.ca /labs/ scl/ spgll/ For M-SBL implementation, we used the original im- 
plementation by David Wipf. Since the mixed norm approach and M-SBL do not provide a exact 
A;-sparse solution, we used the support for the largest k coefficients as a support estimate in 
calculating the perfect recovery ratio. Figs. |2a) and (b) show the recovery rates for = 8 and 256, 
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5 10 15 20 25 30 35 40 



Fig. 4. Sequential CS-MUSIC, Sequential CS-MUSIC without greedy subspace estimation, and CS-MUSIC/SA-MUSIC 
performance for various snapshot numbers. Simulation parameters are n = 128, fc = 8, r = 4 and SNR = 30dB. 

respectively. Sequential CS-MUSIC outperforms S-OMP and the original CS-MUSIC/SA-MUSIC 
consistently, and its performance nearly achieves those of M-SBL and the mixed norm approaches. 
Note that the performance of M-SBL and the mixed norm approaches were identical. Indeed, the 
additional sampling cost for a sequential CS-MUSIC compared to the M-SBL or the mixed norm 
approaches is very small. Considering that any subspace method needs additional redundancy (i.e. 
m > A; + 1) to avoid ambiguity in the signal subspace estimation, we believe that sequential 
CS-MUSIC nearly achieves the optimum performance. Furthermore, this high performance can 
be achieved at negligible computational complexity. Note that the complexity of the sequential 
CS-MUSIC is only a fraction of those of M-SBL and the mixed norm approaches, as shown in 
Figs. I2c)(d) for = 8 and 256, respectively. 

To show the dependency of recovery performance on the condition number of X, we conducted 
simulations for two different types of X. More specifically, the j-th diagonal term of A in Eq. ([T5l) 
is given by aj = t^~^ for j = 1, • • • , r. The results in Fig. [6ta) provide evidence that sequential 
CS-MUSIC is not greatly affected by the condition number of X, and appears less sensitive than 
M-SBL. Next, we performed simulation studies for the different types of RIP conditions using 
various MMV algorithms. More specifically, we assumed that each component of a sensing matrix 
follows J\f{a, 1/m) and then normalized each column of A to have a unit norm. The mean values 
are set to a = and 1, where a larger a represents a worse RIP condition. In this simulation. 
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Fig. 5. Performance of various joint sparse recovery algorithm at (a) A'^ = 8 and (b) A'' — 256, when n — 128, k = 
8, r = 4, and SNR = 30dB. (c)(d) Average CPU time for iV = 8 and = 256, respectively. 



= 64 and the other parameters are the same as before. Fig. Ob) shows that sequential CS- 
MUSIC is more robust that the original CS-MUSIC/SA-MUSIC for unfavorable RIP conditions. 
However, compared to M-SBL, the sequential CS-MUSIC appears less robust to unfavorable RIP 
conditions, which is commonly observed in most of the greedy approaches. 

C. Fourier Sensing Matrix 

Finally, we conducted similar numerical experiments using the Fourier sensing matrix. In this 
case, the source model in Eq. ( fTSl ) is set to be complex valued. Fig. |7] illustrates the recovery 
performance of various MMV algorithms for the Fourier sensing matrix when n = 128, k = 8, 
and r = 4 at the SNR of 30dB for N = 5. We again observed similar performance improvement as 
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5 10 15 20 25 30 35 40 5 10 15 20 25 30 35 40 

m m 



Fig. 6. Recovery rates of M-SBL, sequential CS-MUSIC and CS-MUSIC/SA-MUSIC for (a) different condition numbers 
and (b) various RIP condition of A (larger mean value provides worse RIP condition). The simulation parameters were 
n = 128, it = 8, r = 4, iV = 64 and SNR=30dB. 

in the Gaussian matrix. Note again that the performance of M-SBL and the mixed norm approaches 
were identical. However, compared to the Gaussian cases, a Fourier measurement has redundancies 
in imaginary information, which improves the overall recovery performance. 




Fig. 7. Performance of various joint sparse recovery algorithm from Fourier sensing matrix at when n = 128, k = 
8, r = 4, iV = 5 and SNR = 30dB. 
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VI. DISCUSSION 

A. Relation to Prior Constrained MUSIC and Sequential MUSIC 

In a prior constrained MUSIC algorithm |[T5l . the prior knowledge of a direction-of-arrival 
(DOA) is incorporated into the MUSIC criterion to improve estimation performance by filtering 
out the known source directions via orthogonal projections. If the prior knowledge of a partial 
support is exact, the prior constrained MUSIC is closely related to the generalized MUSIC step 
in CS-MUSIC/SA-MUSIC with an exact partial support estimate. However, as shown in [16| as 
well as in this paper, such an algorithm is affected by the resulting perturbation of the augmented 
subspace if the number of snapshots or SNR is not sufficiently high or the partial support knowledge 
is errorneous. There have been several approaches to improve the noise robustness of a prior 
constrained MUSIC (see |fT6l and references therein); however, to the best of our knowledge, we 
are not aware of any existing method that improves the robustness of an augmented signal subspace 
estimate using sequential subspace estimation and support filtering. 

Sequential MUSIC and its variations 1241, flU, 1261, 1271 may appear closely related to the 
proposed method. Indeed, Davies et al. ll28l showed that the Recursively Applied and Projected 
(RAP)-MUSIC (271 is equivalent to a subspace S-OMP (SS-OMP) step - a partial support recovery 
estimation part of CS-MUSIC. However, as shown in |f7|, the SNR requirement of a subspace S- 
OMP is much tighter than that of a generalized MUSIC step. Therefore, switching from sequential 
MUSIC to the generalized MUSIC step would be beneficial. Moreover, our analysis in this paper 
showed that the sequential subspace estimation further relaxes the SNR condition sequentially. 
This implies that even if the SNR condition of the generalized MUSIC is not satisfied initially, 
during the sequential subspace estimation the SNR condition can be met and the overall recovery 
performance can be improved. To the best of our knowledge, this type of techniques have not been 
reported for any existing sequential MUSIC algorithms IIH, (251, 1261, ll27l . 

B. Sparsity Estimation 

So far, our derivation assumes the prior knowledge of support size. In our previous work Q, 
we derived a sparsity estimation algorithm. Note that this algorithm can be incorporated at each 
greedy step to make the algorithm work, even without knowing the sparsity level a priori. In 
addition, there are various heuristics that could be used to estimate the sparsity in MUSIC type 
parametric methods. However, they have a nonzero probability of being wrong if the measurement 
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is noisy. Typically, these sparsity estimation algorithms tend to overestimate, so various model 
order selection criteria have been often incoiporated to avoid this overestimation [29] . 

Note that the RIP condition (^2fc_^,_|_/ < 1 implies that the maximal sparsity level that our algorithm 
can recover does not exceed m — 1, and the corresponding submatrix Aj for the support estimate 
/ is always full column ranked. This implies that, as long as the condition number of is 
not bad and the noise levels are sufficiently small, we can implement thresholding techniques in 
a reconstruction domain to find a sparse signal, similar to other non-parametric sparse recovery 
approaches like iterative thresholding, M-SBL, etc. When such thresholding scheme may not be 
sufficiently accurate, the current technique has limitations and we need a new way to estimate the 
sparsity level. Though the sparsity estimation is very important topic, this is beyond scope of the 
current work, and will be reported elsewhere. 

C. Limitation of Noisy Analysis 

Even though our noisy analysis provides useful insight on the origin of the noise robustness of 
a sequential CS-MUSIC, current analysis has two limitations. First, the analysis is based on the 
Gaussian sensing matrix using an asymptotic argument. Hence, the analysis should be modified 
for a general sensing matrix such as Fourier. An RIP based analysis in SA-MUSIC |12| would 
work toward this goal. Second, the noisy performance analysis is based on comparing sufficient 
conditions. Since a sufficient condition is often more restricted than necessary, the analysis in this 
paper should be understood as a more conservative comparison. 

VII. CONCLUSION 

In this paper, we derived two greedy strategies to improve the noise robustness of recent hybrid 
joint sparse recovery algorithms such as CS-MUSIC and SA-MUSIC. Although these hybrid algo- 
rithms significantly outperform any other conventional greedy MMV algorithms, the performance 
improvement is reduced for a limited number of snapshots. We showed that the performance 
degradation is due to a perturbation in an augmented signal subspace estimation originating from 
an inaccurate subspace or partial support estimation. Furthermore, we demonstrated that even with 
limited number of snapshots, there are two different ways to improve the noise robustness of 
augmented signal subspace estimation: one by sequential subspace estimation and the other by 
filtering out incorrect support. We further explained that the two greedy steps are byproducts of a 
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novel generalized subspace criterion. Theoretical analysis in noisy situations revealed the origins of 
the noise robustness of the proposed algorithm and led to the identification of sampling conditions 
where each greedy step becomes beneficial. Extensive numerical simulation demonstrated that the 
new algorithm consistently outperforms the existing greedy algorithms and nearly achieves optimal 
performance with minimal computational complexity. 

Appendix A 

To obtain the perturbation bound Eq. ^ in sequential subspace estimation, we use the following 
theorem. 

Theorem 7: ||30l Assume that G € M'"^'? has the singular value decomposition 

G = UT.V* = UiT.iV* + UoTqVq =Gi + Gq 

where d := Ui^iV^ and Gq := UqTqV^ . Also, for a perturbed matrix G E M™^^ of A, assume 
that G has the singular value decomposition 

G = jjtV* = jJitiV^ + f/oSoV? = Gi + Go 

where Gi := C/iSiVi and Gq := UqTqVq , and Ui and Ui (or Uq and Uq) are the matrices of 
same size. If there exist a > and (5 > such that 

o-mm(Gi) >a + 5 and fTmax(Go) < a, 

then for every unitary invariant norm, 

sine(i?(Gi),i?(Gi)) = \\PR^G,)-PRiG,)\\ < ^, 

where 

e:=max(||i?i||,||ii2||), Ri := -{G - G)Vi, R2 := -{G - G)*Ui. 

(Proof of Eq. For a noiseless measurement [Aj^_^^^ S], we have 

fji > • • • > cjfc > C7fc+i = • • • = ak+i = 

so that we have amm{Uit.iV{) > ak{[Ai^_^^, S]) - \\S - S\\ and o-max(f^o5^oVo*) = so that if 

we have A := \\Pr(^s) ~ Pr{s)\\ < '^fc([^/fe_,+, S]), we can apply Theorem|7] Then, we have 

e ^ max{\\{S - S)Vil\\{S - SyUiW) ^ A 

^ B]) - \\S - S\\ - ak{[Aj,_^^, S]) - A 

if A < ak{[AI^_^^^ S]). This concludes the proof. 
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Appendix B 

By the assumption \I n suppX| > k — r + 1 and the generalized MUSIC criterion, we have 
} S])^'?!!^ ~ *-* ^'^^ ™y ^ ^ suppX. Then, for any j ^ suppX and q G suppX, we have 
IIP-L a -IP _ IIP-L a ||2 

> II P-*- a ■ II 2 — II P-*- a II 2 

- IIP-L a -112 _ IIP-L a l|2 

(16) 

since = P for any orthogonal projection operator P, where Ik-r,q C (/ \ {g}) n suppX. For 
j ^ suppX, Aj is statistically independent from R{[Aj\^jj S]) so that n-ll-P^j-j^ 
chi-squared random variable with at least m — k — r + 1 degrees of freedom so that 

liminf min a*^^",,. 5,^a,- > 1 — 7(1 + q), (17) 

where 7 := lim„_j.oo k/m and a := lim„_j.oo r/k. 

On the other hand, for any q G suppX, m||aq|p is a chi-squared random variable of m degrees 
of freedom, so that by Lemma 3 in |[22l we have lim max ||a„|p = 1. Since [Ajj^_^ S] has 

ra— i-oo gGsuppX ''''' 

a full column rank, using the bound in Eq. Q for any Ik-r,q C (suppX n /) \ {q}, we have 

A 

ll^i?([A.,„.„ s]) s])\\ < ^ ^,([^ s])-A- 

Tc(suppXn/)\{g} 

If we let a{q,I) = max{fTfc([A/^ S]) -.T C Ic\ {q}, \T\ = k - r}, where Ic = {qi, ■ ■ ■ 
and 

d-{qi,I) > a{q2,I) >■■■> a{q\j^\,I), 

we have 

A 

II^R([A.,_„,, 5]) - ^iJ([A.,_^„ S])\\ < ^,(5,_^,J)_A ^^^^ 

for g = gi, • • • , ^fc-r- Hence, by ( fT6l ). ( fTTl ) and ( fTSl l. (fT2l) holds if we have 

1 - 7(1 + a) - — ^- > 

ak{qk-r,I) - A 

in the large system limit. This completes the proof. 
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